function [scaled_ATT_bs_itr,ATE_bs_itr,W_bs_Ghat_itr] = generate_delta_CI_bootstrap_cubic(bsimax,sample,wave,tcost,covariate,winter_prevuse,SMC,savedir,fbase)
% This function estimates confidence intervals for the difference between 
% EWM point estimate (from cubic rule) and the scaled ATT.

% Inputs: 
% (1) bsimax: number of bootstrap runs
% (2) sample: cross-sectional dataset merged using elec_wave367_cross.csv
% and propensity scores and renamed as data_estimation_pooled_ps
% (3) wave: indicates whether the output is wave-specific (=3,6,or 7) or
% for the pooled sample (=0)
% (4) tcost: private marginal cost of implementing the program per household
%  >0 represents cost savings; =0 represents kWh reduction 
% (5) covariate: indicates covariate for analysis: income, size, vintage,
% minimum of baseline consumption, maximum of baseline consumption, or
% standard deviation of consumption
% (6) winter_prevuse: indicates whether baseline consumption is calculated 
% as the mean of consumption in winter months (Jan and Feb) or as the mean
% of specified pre-treatment periods.
% (7) SMC: =1 use social marginal cost; =0 use retail electricity price
% (8) savedir: directory to save the savings table 
% (9) fbase: string used to indicate the propensity score and baseline
% months specification for output table 

% Output: 
% (1) Scaled_ATT_bs_itr: array containing scaled att of each bootsrap sample
% (2) ATE_bs_itr: array containing ATE of each bootstrap sample 
% (3) W_bs_Ghat_itr: array containing savings from previously esimtated
% EWM rule on bootstrap sample 
%% Call CPLEX
addpath 'C:\Progra~1\IBM\ILOG\CPLEX_Studio129\cplex\matlab\x64_win64'

% Setup CPLEX optimization options
opt = cplexoptimset('cplex');
opt.parallel = 1;
opt.threads = 4;
opt.simplex.tolerances.feasibility = 1e-08;
opt.mip.tolerances.integrality = 1e-08;
opt.mip.strategy.variableselect = 3;
opt.mip.strategy.nodeselect = 2;
opt.mip.strategy.lbheur = 1;
opt.mip.limits.cutpasses = -1;
opt.display = 'off';

%% Gather inputs
[X,Y,D,n,ps,Yscale,x1_wave,prevuse_wave,Xscale] = generate_inputs_cubic_rule(sample,wave,tcost,covariate,SMC);
opower_paper_wave=sample.opower_paper_wave;

%% Estimate g using IPW
g = Y.*(D./ps - (1-D)./(1-ps)); % elements of W(G) - W(G_0)

%% Solve for rule using cplex 
% first BS permutation- original sample
bsperm = [1:n]';
gn = g(bsperm,:);
Xn = X(bsperm,:);

% number of variables (including the constant)
k = size(X,2);
bsi=0;

tic1= tic;
% Solve for optimal treatment rule
    % generate input for cplex package
    [f,f_n,Aineq_d,Aineq_i,bineq,lb,ub,ctype,Ind] = generate_cplex_inputs(bsperm,g,X,bsi,k,gn,Xn);
    % minimization with treatment decreasing in pre treatment usage
    [sol_pd, v_pd] = cplexmilp(f,Aineq_d,bineq,[],[],[],[],[],lb,ub,ctype,[],opt); %the cplexmilp function is designed for a minimization problem
    % with treatment increasing in pre treatment usage
    [sol_pi, v_pi] = cplexmilp(f,Aineq_i,bineq,[],[],[],[],[],lb,ub,ctype,[],opt);
    
    % get beta and v 
    if (v_pd < v_pi)
            beta = sol_pd(1:k,:);
            v    = v_pd;
        else
            beta = sol_pi(1:k,:);
            v    = v_pi;
    end
    % get solved assignment in_Ghat and cost_saving
    in_Ghat=(X*beta>0);
    this_cost_saving=nanmean(g.*in_Ghat)*Yscale;

fprintf('cplex solution takes %.2f sec\n',toc(tic1))

%% Bootstrap for delta_CI

% initiate outputs
scaled_ATT_bs_itr = zeros(bsimax,1);
ATE_bs_itr = zeros(bsimax,1);
W_bs_Ghat_itr = zeros(bsimax,1);

% specification suffix 
if (tcost)
    if (SMC)
        s_speci = 'smc';
    else
        s_speci = 'pmc';
    end
else
    s_speci='kwh';
end
% shared name for each bootstrap batch
s_now_batch=datestr(now(),'yyyymmdd_HHMMSS');
% Create a saving directory for the i_BS files
 fpath_BS = ['.\intermediate_data\EWM_results\BS\cubic\',...
     s_now_batch,'_',covariate,'_',s_speci,'_','deltaCI','\'];
if ~exist(fpath_BS,'dir')
   mkdir(fpath_BS); 
end

sprintf('bootstrap:%s%s%d',s_now_batch,covariate,tcost)

for bsi = 1:bsimax 
    tic1=tic;
    rng(bsi);
    dt_start = now();
    % DRAW BOOTSTRAP SAMPLE WITH INDEXES bsperm
    bsperm = randi(n,[n 1]); % DRAW BOOTSTRAP INDEXES FOR THE NEXT SAMPLE
    
    %% obtain att, ate, and savings from applying the original G_hat on bootstrap sample
    [scaled_ATT_bs,ATE_bs,W_bs_Ghat] = ewm_rct_test_cubic_input(D,Y,ps,beta,X,in_Ghat,g,Yscale,bsperm,n);

    scaled_ATT_bs_itr(bsi,1) = scaled_ATT_bs;
    ATE_bs_itr(bsi,1) = ATE_bs;
    W_bs_Ghat_itr(bsi,1) = W_bs_Ghat;
    
    %% save each BS result
    dt_outp = now();
    dt=[dt_start,dt_outp];
    save_BS_delta_CI([fpath_BS,sprintf('%s_BS_%d.mat',s_now_batch,bsi)],...
        bsi,s_now_batch,in_Ghat,scaled_ATT_bs,ATE_bs,W_bs_Ghat,bsperm,dt);
    % progress markers
    fprintf('Boot #%d took %2.3f sec\n',bsi,toc(tic1))

end


end